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ABSTRACT 

This paper considers the effects of turbulence on mean motion resonances in ex- 
trasolar planetary systems and predicts that systems rarely survive in a resonant con- 
figuration. A growing number of systems are reported to be in resonance, which is 
thought to arise from the planet migration process. If planets are brought together and 
moved inward through torques produced by circumstellar disks, then disk turbulence 
can act to prevent planets from staying in a resonant configuration. This paper studies 
this process through numerical simulations and via analytic model equations, where 
both approaches include stochastic forcing terms due to turbulence. We explore how 
the amplitude and forcing time intervals of the turbulence affect the maintenance of 
mean motion resonances. If turbulence is common in circumstellar disks during the 
epoch of planet migration, with the amplitudes indicated by current MHD simulations, 
then planetary systems that remain deep in mean motion resonance are predicted to 
be rare. More specifically, the fraction of resonant systems that survive over a typical 
disk lifetime of ~1 Myr is of order 0.01. If mean motion resonances are found to be 
common, their existence would place tight constraints on the amplitude and duty cycle 
of turbulent fluctuations in circumstellar disks. These results can be combined by ex- 
pressing the expected fraction of surviving resonant systems in the approximate form 
'Abound ~ C'/A'orb^^^, where the dimensionless parameter C ~ 10 — 50 and A'orb is the 
number of orbits for which turbulence is active. 

Subject headings: MHD — planetary systems — planetary systems: formation — plan- 
ets and satellites: formation — turbulence 



1. INTRODUCTION 



An increasing number of the observed extrasolar planetary systems have been discovered to 
contain multiple planets. A growing subset of these multiple planet systems have period ratios 
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close to the ratios of small integers and hence are candidates for being in mean motion resonance 
(e.g., Mayor et al. 2001, Marcy ct al. 2002, Butler et al. 2006). In true resonant configurations, the 
orbital frequencies not only display integer ratios, but, in addition, the relevant resonant arguments 
(angular variables composed of the osculating orbital elements) display constrained oscillatory time 
evolution (for further discussion, see Murray & Dermott 1999; hereafter MD99; see also Peale 1976). 
As a result, these resonant configurations represent rather special dynamical states of the planetary 
system and their existence (if/when observed) places interesting constraints on the formation and 
long term evolution of these systems. More specifically, as we show in this paper, planetary systems 
are easily knocked out of mean motion resonance. One source of perturbations acting on the planets 
is turbulent fluctuations in the circumstellar disks that originally formed the planets. Since these 
disks are also necessary ingredients in the current picture of planetary migration (e.g., Papaloizou 
k. Terquem 2006), and since coupled migration seems necessary to put planets into mean motion 
resonance (e.g., Lee & Peale 2002; Lee 2004; Moorhead & Adams 2005; Beauge et al. 2006), the 
effects of disk turbulence on the resonances is important to understand. 

This paper analyzes the effects of turbulence on planets in or near mean motion resonance 
using both a semi-analytic treatment (§2) and direct numerical integrations (§3). The semi-analytic 
approach makes extreme approximations in order to obtain tractable equations, but it allows for an 
explicit determination of how the long term behavior depends on the basic physical variables (e.g., 
planet masses, semimajor axes, type of resonance, turbulent forcing amplitude, etc.). In contrast, 
numerical integrations can be carried out to high accuracy including all 18 dynamical variables of the 
problem (for 3-body systems). Although the simulations demonstrate what the long term behavior 
will be, they don't automatically specify the relationships between the aforementioned physical 
variables. Fortunately, both approaches give consistent results for the intermediate time scales of 
interest. In particular, we find that the expectation value of the effective energy of the resonance 
- considered as a nonlinear oscillator - increases linearly with time. Similarly, the probability that 
systems stay in resonance decreases with time. For the stochastic pendulum model, the fraction 
of systems in resonance decreases as the square root of time for the simplest case in which the 
systems freely random walk in and out of a bound state; the decay is exponential if systems cannot 
return to resonance after leaving. The full numerical integrations give results that are intermediate 
between these extremes. For the model equations, the leading coefficients in the time-dependent 
relationships are given by the amplitude of the fluctuations, the time intervals of their forcing, and 
the natural oscillation frequency of the libration angle of the original resonance. This paper thus 
provides a relatively simple description of the effects of turbulent fluctuations on planetary systems 
in or near mean motion resonance. 
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2. SEMI-ANALYTIC MODEL EQUATIONS 

2.1. Pendulum Model for the Resonance 

For the sake of definiteness, we consider the case in which a larger planet on an external orbit 
perturbs a smaller planet on an interior orbit. The larger body is then assumed to have orbital 
elements that do not change with time (or at least vary much less than those of the smaller inner 
planet). We are thus making the most extreme approximation and keeping only the leading order 
terms in order to obtain an analytic description. Following MD99, the equation of motion for the 
libration angle 93 for two planets in resonance reduces to that of a pendulum, i.e., 

+ WqSiik/? = 0, (1) 

where the natural oscillation frequency ujq is given by 

ujI = -3j|CrJ^el^4l . (2) 

Here, Vt and e are the mean motion and eccentricity of the inner planet. The integers j2 and 
depend on the type of resonance being considered. The parameter C,., which is taken to be constant, 
is given by 

Cr = fi^afdia) , (3) 

where fj, = rUo^t/M^ and where mout is the mass of the outer (perturbing) planet. The quantity 
afd{a) results from the expansion of the disturbing function (MD99), where the parameter a is the 
ratio of the semimajor axes of the two planets, i.e., a = 01/02. This approximation neglects terms 
of order 0{fj,). Finally, we note that more complicated analytic models for mean motion resonances 
can be derived (e.g., Holman & Murray 1996; Quillen 2006), but they are qualitatively similar to 
the pendulum equation considered here. 

For much of this analysis, we are interested in the 2:1 mean motion resonance since observed 
extrasolar planetary systems are (sometimes) found in or near such a configuration, and since these 
resonances are generally the strongest. For this case, the integers j2 = —1 and j4 = —1. For the 2:1 
resonance, afdia) ~ —3/4, where a is the ratio of semimajor axes of the two planets. With these 
specifications, the natural oscillation frequency uq of the libration angle is given approximately by 

Lol^^fiQ'e. (4) 

In typical cases, where planet masses are 1000 times smaller than the stellar masses, and the 
eccentricity e ~ 0.1, we find that loq/Q ~ 10~^. As a result, the period of oscillation for the 
libration angle is typically ^ 100 orbits. 

For completeness, we note that for first order resonances, the oscillation frequency often 
has additional contributing terms (MD99). The order of the additional terms differs from that of 
equation ([2]) by a factor of order 0(/i/e^), where e is the eccentricity of the inner planet. Under 
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some circumstances of interest, /i ~ 10~ and e ~ 0.1, so this correction can be of order unity. We 
are thus considering only the simplest form of the pendulum equation in this analysis, in order to 
gain a general understanding of the effects of turbulence on resonant systems. In principle, however, 
care must be taken to include all of the relevant terms. 



2.2. Turbulent Perturbations 

Given the pendulum model for the resonance (from MD99), the next step is to include the 
effects of additional perturbations produced by turbulence in the circumstellar disk material that 
remains after the planets have been formed. Since this disk is generally thought to induce planetary 
migration, some disk material is likely to be present, for a typical time scale of 1 - 10 Myr. Under 
many circumstances, these disks are expected to be turbulent (Balbus & Hawley 1991) due to the 
magneto-rotational instability (MRI). Although MRI can be shut down if the ionization rate is 
not high enough (Gammie 1996), star forming regions can experience enhanced cosmic ray fluxes, 
and hence enhanced ionization rates, due to cosmic rays from supernovae becoming trapped in 
the magnetic fields of molecular clouds (Fatuzzo et al. 2006). In this work, we consider turbulent 
fluctuations to be present and consider our previous model (Laughlin et al. 2004; hereafter LSA) 
to provide the parameters of the perturbations. Note that the LSA model describes turbulence 
is disks with no gaps, applicable to low-mass embedded planets, whereas this paper primarily 
considers giant planets that will clear gaps. As a result, we must modify the LSA formalism to 
include this effect (see below). 

In principle, the turbulent fluctuations in the disk can affect either planet, i.e., the outer 
perturbing planet (in this analysis) or the inner perturbed planet. If the outer planet is subject 
to turbulent fluctuations, they have the effect of shaking the effective base of the pendulum. As a 
result, they modify the effective potential energy of the pendulum, and the stochastic term enters 
into the equation of motion as follows: 

^ + K' + e(t)]sin(/. = 0, (5) 

where ^(t) is the stochastic forcing term. This is the form expected when the circumstellar disk 
surrounds both planets and primarily influences the outer one. In contrast, if the disk (the torques 
produced by turbulence in the disk) primarily influences the inner planet, then the fluctuations 
provide an additional acceleration term and the equation of motion takes the form 

+'^oSin(/? = -^(t) . (6) 

For the sake of simplicity, we will primarily consider the first case (eq. [5]). Notice also that both 
types of fluctuations can be included. 

The inclusion of turbulent fluctuations thus provides a stochastic forming term in the pendulum 
equation of motion ([T|). The forcing term ^(t) can be considered as a series of "kicks" produced 
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by short-term forcing due to turbulence. To completely specify this forcing effect, one needs to 

determine the spectrum of fluctuation amplitudes, the timing of the fluctuations, and any possible 
correlations in the fluctuation timing. This latter property determines the so-called "color" of the 
noise (fluctuations) and can have (perhaps) surprising effects on the long-term evolution of the 
stochastic pendulum and hence the resonance angles of the planetary system (Mallick k, Marcq 
2004, hereafter MM04). 

The formulation for turbulence developed in LSA shows that the turbulent fluctuations are 
reloaded into the disk on roughly an orbital time scale. In order to consider the turbulent torques 
to be fully independent, we must use the longest relevant orbital period, which corresponds to that 
of the outer disk edge (in the calculations of LSA). This period at the outer disk edge is about 
twice that of the planet being forced. Since the outer (perturbing) planet is being forced by the 
turbulence in the present context, and since the outer planet is in a 2:1 mean motion resonance 
with the inner planet, this turbulent forcing period is approximately four times the period of the 
inner planet. 

Next we need to account for the fluctuations in the equation of motion. To the same order of 
approximation used to derive the pendulum equation for the resonance, the time rate of change of 
the mean motion is given approximately by 



where J is the angular momentum of the inner planetary orbit. Here we consider the turbulent 
fluctuations to provide discrete changes in the angular momentum on time scales comparable to 
the orbital period Pd at the disk edge, where (as discussed above) we expect Pd ~ ^Pm-h = Stt/Q,. 
Thus, the time scale for independent stochastic perturbations to act is given by At ~ Svr/Q. 
Since cPip/dt'^ ~ dfl/dt (see equation [8.39] from the derivation of MD99), the forcing term .^^ oc 
Q,{J/J) ~ Q,[{AJ)k/ J]6{[t] — At). As a result, we can write the stochastic differential equation in 
the form 



In this equation, [(AJ)/J]fe represents the fractional change in the angular momentum of the planet 
for a given realization of the turbulence. The subscript 'A;' indicates that this increment of angular 
momentum is delivered at discrete time intervals, which can be counted and labeled with an integer. 
To account for this timing, the Dirac delta function is periodic with period At, so that the quantity 
[t] is the time measured mod-At. 

The amplitude [(AJ)/J]k is calculated in LSA for disks without gaps (see also Nelson & 
Papaloizou 2003, hereafter NP03, and the discussion below). The resulting torques produce angular 
momentum perturbations with zero mean and well-defined RMS amplitudes of order [(AJ)/J]rms = 
jrms ~ 0.004. Other three dimensional MHD simulations produce fluctuations with comparable 
amplitudes (e.g.. Nelson 2005). 




(7) 




-6- 



For planets that clear gaps in the disk, as considered here, we must include a reduction factor 
Tji to account for the absence of the full complement of disk material near the planet. In practice, 
we find that Tf{ ~ 0.1, where this value can be estimated in multiple ways. First, we note that 
numerical simulations of gap formation in turbulent disks (NP03) show that the gaps are not 
completely cleared out; instead, the surface density at the gap minimum is reduced by a factor 
of ~ 20 from the unperturbed level (see Fig. 1 of NP03). Thus we expect the reduction factor 
to lie in the range 0.05 < F/j < 1. Next we note that material from both the outer part of the 
gap (where the reduction in surface density is less severe) and outside the gap (where there is no 
reduction) provide some contribution to the torques. Indeed, the dominant contribution of torques 
is expected to come from intermediate length scales (see Johnson et al. 2006), specifically scales 
larger than the disk scale height, but not too much larger (given that gravitational forces decrease 
with distance). To model this effect, we start with the formalism of LSA and remove disk material 
corresponding to a quadratic gap profile with zero surface density at the center (e.g., Goldrcich & 
Sari 2003) and then calculate the reduction factor F^j as a function of gap width. Let Wg denote the 
half-gap width, i.e., the location on cither side of the planet where the quadratic gap profile joins 
onto the unperturbed surface density distribution. The planet is assumed have semi-major axis 
Up, which also defines the location of the gap center, where the surface density is taken to vanish. 
Then we find that Tr 0.2, 0.1, and 0.05 for total gap width 2wg/ap = 1, 1.5, and 2, respectively. 
Although the exact value of the reduction factor will depend on the details of the gap shape, we 
use Fff = 0.1 as a representative value for this work. We also note that these systems have two 
planets, so that each planet can clear disk material, but each planet can also experience torques; 
we assume here that these two competing effects cancel. 

For completeness, we note that the surface density structure in the local neighborhood of a 
massive planet can be strongly perturbed, with well-defined wakes and even a circumplanetary disk 
within the gap (e.g., see Fig. 10 of NP03). Although global transport seems to be unaffected by 
these distortions, their effects on stochastic torques remain unclear and should be considered in 
future work. 

Since the forcing amplitude [{AJ)/J]k plays an important role, it is useful to have an heuristic 
understanding of the expected value. For a circunistellar disk with surface density S, the torque 
exerted on a planet can be written in terms of the benchmark value To = 27rGSrmp. The torque 
due to turbulence is expected to be a fraction Jt of this physical scale, where numerical simulations 
suggest that /t ~ 0.05 (e.g., Johnson et al. 2006; Nelson 2005). This torque acts over a time 
interval to produce a change in angular momentum. As argued above, this calculation uses a 
time interval of four orbital periods, long enough so that the turbulent torques can be taken to be 
independent. The torque thus acts on a timescale 4Por6 = 87r(r^/GM*)^/^ and the net change in 
angular momentum (A J) = 4:PorbfTTD- Next, we include the reduction factor F^ due to the disk 
gap (where we expect F^j ~ 0.1). Finally, the orbital angular momentum J of the planet is given 
by J = mp[GM,r(l - e2)]V2. We can absorb the eccentricity dependence into the definition of fx- 
Putting these results together, we can write the fractional change in angular momentum over one 
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time interval (4Porb) in the form 

where the second equahty is scaled for a 1 AU orbit and a typical surface density at that radius. 
Note that we expect the surface density profile to have a nearly power-law form S r^^, where 
p ~ 3/2, so the quantity Sr^ increases slowly with radius. Since giant planets form at larger radial 
location and migrate inwards, the torques will be somewhat larger during the earlier phases of 
migration. 



2.3. Time Evolution of the Stochastic Pendulum 

With the mean motion resonance modeled by a pendulum equation, and with the character- 
istics of the turbulent fluctuations specified, we now find the time evolution of the system. If we 
define a dimensionless time variable 

T = UJot, (10) 

where loq is the oscillation frequency of the libration angle, and is defined by equation ([4]), the 
equation of motion takes the form 

^ + [l + rjkH[r]-^r)]sm^ = 0, (11) 

where 6{x) is the Dirac delta function. In this formulation, At is the time interval between 
stochastic kicks, the dimensionless time variable is measured mod-Ar, and the amplitude of the 
kicks is given by 

1 n /AJ\ _ 2 /AJ\ 

The quantity r/^ is thus a stochastic variable that has a distribution of values inherited from the 
distribution of values of [(AJ)/J]fc, which in turn is determined by the properties of the turbulence. 
The kicks in angular momentum can be either positive or negative, with no expected asymmetry, 
so that (r/fc) = 0. We can thus characterize the amplitude by the RMS (root-mean-square) of the 
variable. Using the amplitudes expected from numerical MHD simulations (§2.2), this amplitude 
is expected to lie in the range 

r?rms = (r?I>'/'~ 0.002 - 0.009, (13) 

where we have used typical values of the mass ratio fj, = 10^^ and the eccentricity e = 0.1 to eval- 
uate the result. The typical expected value of the RMS fluctuation amplitude is thus 77rms ~ 

0.005. 

These acceleration kicks occur at typical time intervals At ~ Svr/O, which sets the corresponding 
dimensionless time interval to be At ~ Sirujo/i}. In general, both the amplitudes and the time 
intervals between forcing kicks will vary from cycle to cycle according to their (respective) distri- 
butions. However, we can scale out the time variations (Appendix A; Adams & Bloch 2008) and 
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characterize the turbulence with a well-defined (single-valued) forcing time. As a result, we present 
the remainder of this analysis using a single forcing time interval and focus on the effects of the 
distribution of forcing amplitudes r]k- 

To analyze the behavior of the stochastic pendulum, as defined above, it is Tiscful to work in 
terms of phase space variables. We thus rewrite the equation of motion into two parts, i.e., 

^ = V, and ^ = - [1 -|-?7l sinw. (14) 
ar ar 

We then define the probability distribution function for the phase space variables, P{(p, V; t), which 
obeys the Fokker-Planck equation (MM04; Binney & Tremaine 1987) 

dP ^BP . dP V . 2 d^P 

_ + y__smc,— = -sm (15) 

where the phase space diffusion constant T> is set by the amplitude of the fluctuation spectrum. 
Specifically, we can write the diffusion constant in terms of previously defined variables, 

P^M)~M)^1 (16) 

At Sttloq ' ^ ' 

Next we argue that the libration angle (p varies rapidly compared to the velocity V. This claim can 
be verified by noting that dV/dr ~ r/fe so that V ~ r^/^ (since the random variable rjk implies a 
random walk growth). However, the libration angle obeys the equation d^p/dr = V, which in turn 
implies that cp ~ t^^"^ (see MM04 for further discussion). In the long term, the libration angle thus 
varies faster than V and we can average the Fokker-Planck equation over the angle (p to obtain the 
time evolution equation for the averaged probability distribution function p{V; r), i.e.. 



This equation has the well-known solution 



Dt 



(18) 



where we have defined an effective diffusion constant D = 2X'(sin^(^). In the long time limit, 
(p ^ r^/^, and the angle fully samples the range [0, 27r], so that {sin^ (p) = 1/2 and hence V = D. 
Notice also that the velocity V grows as F ~ t^^'^ in the long time limit. Since the energy of the 
pendulum grows large at long times, the kinetic term dominates the potential energy term, and 
the energy E « V"^ /2. With this substitution, the probability distribution function for the energy 
takes the simple form 

2E' 



P{E;t) = 1 
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(19) 



This result implies that the expectation value of the energy grows linearly with time in the long 
time limit, i.e., 

{E) = ^T. (20) 
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Fig. 1. — Mean time evolution of the energy of a stochastic pendulum with noise terms of various 
"colors". For each curve/line, 10^ realizations of the time evolution of the oscillator have been 
averaged together. The curves show the results for white noise (solid), pink noise (dotted), and 
brown noise (dashed). The first two types of noise lead to energy evolution that is linear in time, 
[E) ~ t, whereas brown noise leads to energy that increases {E) ~ \/t. 
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This linear time dependence can be understood in terms of a random walk of velocity, as shown in 
Appendix B. 

This time behavior is seen in numerical simulations of the stochastic pendulum equation. 
Figure [1] shows the time development of the mean energy, averaged over 10^ realizations, for three 
different types of fluctuations. The solid curve shows white noise, where the fluctuation amplitude 
and the time interval At between kicks are uniformly distributed about well-defined mean values. 
The dotted curve shows the case of "pink" noise, in which the distribution of time intervals is taken 
to have a decaying exponential form; this distribution leads to more time intervals near the lower 
end of the spectrum. As long as the mean time interval and mean fluctuation amplitude remain 
the same, and as long as both are uncorrelated with previous samples, the linear time dependence 
predicted by equation (|20p holds. For completeness. Figure [T] also shows the case of "brown" noise, 
in which the time intervals undergo a random walk (akin to brownian motion) so that a correlation 
of time intervals is present; in this case, the time development of the energy is slower with (E) oc t^^^ 
(see also MM04). Thus, the long term evolution of an ensemble of stochastic pendulums obeys the 
analytic expectations derived above, provided that (1) the fluctuations are independent, (2) the 
long-time limit has been reached, and (3) the system can freely random walk back into a resonant 
(bound) state (from a higher energy, unbound state). 



2.4. Upper Limit on the Fraction of Bound States 

Given the time evolution of the distribution function, we now consider the fraction of low 
energy states as a function of time. Here, states of the system with sufficiently low energy are 
bound and the pendulum oscillates — instead of circulates — so that the system is in resonance. 
We are thus finding an estimate of the fraction of systems that remain in resonance as a function of 
time, where systems leave resonance due to exposure to turbulent forcing. Note that our treatment 
thus far provides only an upper limit to the fraction of bound states because we allow the system 
to freely random walk both in and out of resonance. The fraction of bound states derived in this 
section thus represents the largest possible value. In the following section, we find a more accurate 
estimate by accounting for the reduced probability of systems returning to resonance after leaving. 

As formulated here, the energy E is the kinetic energy, so that the criterion for the resonance 
to be undone (for the libration angle to circulate) is given (approximately) hj E > 1. Note that the 
system must have > 2 to be freely circulating, but systems with energies in the range 1 < E < 2 
have such large libration angles as to not be in a well-defined resonant state. As encapsulated 
in the probability distribution function for the energy, the system always has some probability 
of remaining in resonance, even as the mean energy (E) grows ever larger. This probability is a 
decreasing function of time and is given by 
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e'' dz = Erf (zo) , (21) 
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Fig. 2. — Time evolution of the probability for the planetary system to stay in a mean motion 
resonance as a function of time. The solid curve shows the result of 10"^ realizations of the stochastic 
pendulum equation; in this set of numerical experiments, the system is allowed to come in and out 
of resonance freely. The dashed curve shows the expected asymptotic behavior of the system, i.e., a 
survival probability of the form Abound oc t^^^^. The dot-dashed curve shows the survival probability 
for when a one-way barrier is imposed so that systems that leave resonance (the pendulum becomes 
unbound) are not allowed to re-enter the bound configuration. The dotted curve shows the analytic 
approximation to this probability evolution derived in the text. 
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where zq = . Here Erf(z) is the error function (e.g., Abramowitz & Stegun 1970), which 

can be expanded in the hmit of small z as follows: 

Erf (z) = ^(z-- + — - — + ..) , (22) 
^ ^ ^^F V 3 10 42 / ' ^ ^ 



where small z values correspond to late times. To leading order, the probability that the planetary 
system remains in resonance is thus given by the relation 

^'°"'^'"(^) ^(^) ^Vrjntf ^V1^!VW^' ^'^^ 

where this expression is valid only for sufficiently late times when r > 2/D. For planetary systems, 
the natural clock is set by the orbit time, so in the final equality we measure time using the number 
of orbits iVorb = f^t/(27r). If we then insert the expected fluctuation amplitude T^rms 

^ 0.005 (eq. 

[13]), the above expression simplifies to the form Pbound ~ 640/A^orb^'^^j where A'orb includes only 
the orbits (time) for which turbulence is active. For typical systems, where the number of orbits 
-^orb = 10® — 10®, this expression implies that the probability of remaining in resonance lies in the 
range "Pbound ~ 0.06 — 0.60. Keep in mind, however, that this derivation assumes that planetary 
systems can freely random walk both in and out of resonance, so that equation (1231) represents an 
upper limit to the fraction of surviving resonant systems. But even this upper limit demonstrates 
that turbulence significantly reduces the expected number of resonant systems. In the more realistic 
case, where systems can random walk out resonance more easily than they can random walk back 
into resonance (§2.5 and §3), the fraction of surviving systems will be much lower. 



Equation (j23p shows that the probability of a planetary system remaining in resonance should 
decrease as the square root of time. This dependence (found analytically) is observed in numerical 
simulations using the stochastic pendulum equation. Figure [2] shows how the probability of staying 
in resonance decreases with time for an ensemble of stochastically driven oscillators. The solid 
curves show the percentage of systems remaining in resonance as a function of time, using a sample 
of 10^ systems. The smooth dashed and dotted curves show the analytic predictions for the time 
dependence derived above. 

This estimate of the probability for remaining bound is approximate in several respects. First, 
we have not included the potential energy term [loq cos ip) in the calculation of the probability 
distribution P{E;t). This approximation leads to an uncertainty of ~ \/2 in the estimates pre- 
sented here. More importantly, we again stress that the analysis presented thus far represents an 
upper limit on the fraction of systems that can remain in resonance. In the stochastic pendulum 
formalism, we are implicitly assuming that the system can random walk both in and out of a bound 
configuration, i.e., that the planetary system can random walk in and out of resonance. The model 
equation developed here only has one degree of freedom ((^), so this behavior is natural. In actual 
planetary systems, however, other physical variables are relevant and must have the proper values 
to allow the planets to be in resonance. This difficulty implies that once a planetary systems leaves 
resonance, it will be unlikely to random walk back into resonance. This behavior, in turn, implies 
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that the fraction of planetary systems remaining in resonance will fall below the fiducial law derived 
here. We consider this complication in the following sections. 



The discussion presented thus far assumes that the stochastic oscillator can diffuse in and out 
of a bound state, i.e., into resonance as well as out of resonance. For the simple model case, which 
has only one variable, this behavior is sensible. For real planetary systems, however, the conditions 
of resonance are more complicated than the simple pendulum model. In particular, for highly 
interactive systems (e.g., the observed 2:1 resonance in GJ876; see §3), many orbital elements of 
both planets must have the proper values for the planets to be in a mean motion resonance. In 
such systems, while it remains possible for fluctuations to compromise the resonance, it will be 
unlikely for the system to random walk back into resonance. When a real planetary system is 
knocked out of libration and into circulation, then the chances of returning are reduced because 
the chances for close encounters that drastically change the orbital parameters are significantly 
increased. The fraction of systems that remains bound will thus decrease with time more sharply 
than suggested by the considerations of the previous section. As a result, the model derived above, 
where T^bound oc t~^/^, represents an upper limit on the fraction of systems that remain in resonance. 
In this section, we explore a more complicated model for the probability evolution that includes the 
one-way nature of this process. We also consider the intermediate case of a partially transmitting 
barrier, where systems can return to resonance after leaving, but with reduced probability. 



We now consider a simple generalization of the time evolution of the probability. At any given 
time, the fraction of systems that would remain bound in the absence of the one-way barrier is 
given by "Pbound from equation (f23l) . We then assume that some fraction of these systems are 
"lost" at a well-defined rate, so that the time evolution of the fraction of bound states (in the 
absence of the diffusion term) would be an exponential decay. This ansatz thus assumes that the 
remaining systems are distributed across the possible bound energy values, and that each system 
has some probability of leaving its bound state per unit time. With this assumption, the Fokker- 
Planck equation for the averaged probability distribution function p{V; r) now takes the modified 
(heuristic) form 



where we assume that the bound fraction has the asymptotic form derived above and where the 
parameter 7 encapsulates the details of this "decay process". Note that a similar form arises 
(Caughey 1963) when dissipation is included in the Fokker-Planck equation (these connections 



2.5. Evolution with a One-way Barrier 



2.5.1. Heuristic Model Equation 




(24) 
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should be explored further). Equation (j24|) can be solved to obtain the result: 

= 737^377? exp -— exp "7 ( ^ 



(25) 



With this complication, the fraction of systems that remain bound as a function of time now takes 
the form 

Abound - —FT <3Xp 



ttDt 



-7 



/32r 



(26) 



Although heuristic, this model captures the basic time behavior for stochastic pendulums with a 
one-way barrier, as shown in Figure [2j 



2.5.2. Solutions from Separation of Variables 



For the simplest case of a constant diffusion constant, the solution for a one-way barrier can be 
found using separation of variables. In this case, the solution to the diffusion equation is assumed 
to have the general form p{t, V) = T{t)g{V); the diffusion equation then can be written as 

IdT A 1 <fg 



A' 



T dt D g dV^ 
where is the separation constant. 
The solutions g{V) have the form 

gn(y) = AnCOsAnV , 



(27) 



(28) 



where the subscript denotes one of a series of solutions that satisfy the boundary condition. We 
we invoke the boundary condition that g{V = \/2) = 0, which implies that the distribution func- 



tion vanishes when the kinetic energy reaches E = 
eigenvalues A„, i.e., 



An = ^(n + l/2), 

where n is an integer. The full solution thus takes the form 

oo 

p{t, V) = '^An cos(Any) exp 

n=0 



1. This condition thus specifies the 

(29) 



(30) 



For the initial condition p{t = 0,V) = S{V), where 6{x) is the Dirac delta function, the constants 
are specified so that 

An = 1/V2. (31) 
The fraction of bound states at any time is then given by 



nound = 2 / pit, V)dV = V 
■^0 n=0 



2(-l)" 
7r(n+l/2) 



exp 



--7r2(n + l/2)2t 



(32) 
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This solution, though approximate, indicates that the long term evolution of the ensemble of 
resonant states follows an exponential decay. For the simple solution of equation ([30|) , each "mode" 
decays exponentially, albeit with a differing decay rate. In the long run, however, only the lowest 
order mode survives, and the time evolution becomes purely exponential. 

2.5.3. Further Complications 

The evolution of the actual stochastic pendulum differs from that described by a simple dif- 
fusion equation in several respects. Most importantly, the diffusion constant D is not really a 
constant. In the formulation developed here, D oc sin^ 99, and we have taken the average value 
(sin^ Lf) = 1/2. However, this limit only applies in the long time limit when most of the oscillators 
in the ensemble have large energy, so that their angles uniformly sample the possible values. For 
bound states, by definition, the energies are low, and hence this limiting form is not strictly applica- 
ble. In particular, the typical amplitudes of the bound oscillators will be smaller than average, and 
the true effective diffusion constant will be less than that of the long time limit. In addition to its 
lower value, the diffusion constant will also have some type of time dependence: In the beginning, 
when all of the oscillators in the ensemble have small amplitudes, D will be correspondingly small. 
As time increases, even the bound oscillators will have larger amplitudes and hence the evolution 
is properly described by larger values of D (but still smaller than that of the long time limit). 

2.5.4. Partial Barriers 

As shown above, the time evolution of the number of bound states decays exponentially for 
systems with a strict one-way barrier, i.e., for the case in which systems that leave resonance are 
not allowed to return. We now consider the case in which a given system has some probability 
Pret of re-entering a bound (resonant) state after leaving, where < pict ^1- In this case, all of 
the systems are allowed to stay in play, but only a fraction of the systems can make the transition 
from an unbound state to a bound state. As shown below, even when the fraction prct ^ 1> this 
modification makes a large qualitative change in the long-term behavior of the system. 

We first note that at late times, essentially all of the bound states in the original problem 
(without a barrier) are those that have returned to resonance after leaving. This claim follows 
directly from the above results: The fraction of bound states of all types decreases as ~ 
while the fraction of bound states that have never left resonance decreases as ~ exp[— 7t], where 
the decay rate 7 depends on the details of the diffusion process. At late times, however, the power- 
law behavior wins, and hence most bound states are due to systems that have returned from higher 
energy states. If the higher energy states are allowed to re-enter resonance with probability pret; 
then the fraction of systems in a resonant state will again follow the law, but with a smaller 
normalization. Specifically, the normalization is reduced by the factor p^et and the long-term 
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solution for the fraction of bound states becomes 

/ 8 

Abound ~ Pret ( j • (33) 

This behavior is illustrated in Figure O which shows the time evolution for an ensemble of 
stochastic pendulums. The upper two curves show the evolution for the case in which the systems 
can freely re-enter a resonant (bound) state after leaving. The bottom two curves show the fraction 
of bound states for the case in which the systems are allowed to have only a 20% probability of 
returning to resonance after leaving. The solid and dashed curves show the results of the actual 
numerical integrations, whereas the dotted curves show the analytic approximations derived above. 
Notice the good agreement between the predicted asymptotic forms and the numerical results. The 
fluctuations of the numerically determined fractions about the analytic predictions have amplitudes 
that are consistent with root-A^ noise; note that the numerical ensemble contains only systems 
and the number of bound states is much smaller by the end of the time interval shown in Figure [3j 

On a related note, one can consider the problem in which systems leaving a bound state have 
some probability px of never returning to resonance. Perhaps surprisingly, the long term behavior 
for this case is wildly different than for the previous case. Here, for any nonzero value of px, the 
time dependence of the number of bound states is always an exponential decay. 



2.6. Other Resonances 

As outlined above, the 2:1 resonance is generally the strongest, and planets in such a config- 
uration should survive longest in the presence of turbulence. For other resonances, the oscillation 
frequencies are generally smaller so that the periods of libration are longer. For the case of the 3:1 
resonance, for example, ones finds (MD99) that 

Lol ^ -OmSGnn'^e^ . (34) 

The negative sign indicates that the stable equilibrium point of the oscillator occurs at 99 = vr 
rather than if = 0. After defining (p = (p — ir, the equation of motion for (p becomes the same as 
before. The above frequency is thus lower than that of the (simplified) 2:1 resonance by a factor of 
~ (3/e)^/^. The potential energy of the oscillator is thus less deep by a factor (3/e)^''^, about 5-6 
for typical cases, so that the system is more easily bounced out of a resonant state by turbulence. 

In order to assess the probability of remaining in a bound (resonant) configuration, one can 
evaluate the equations derived in §2.3 - 2.5 using the proper frequency given by equation (j34p . 
For example, the resulting estimate for the probability of remaining in resonance has the form of 
equation (j23p . with the leading coefficient smaller by a factor of (3/e)^/^ ~ 5. In other words, at 
any given time, the probability of a planetary system remaining in a mean motion resonance, in 
the face of turbulence, is inversely proportional to the relevant value of loq. For typical values of 
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Fig. 3. — Time evolution of the fraction of planetary systems that stay in a mean motion resonance 
as a function of time. The solid and dashed curves show the result of 10^ realizations of the 
stochastic pendulum equation. For the top (solid) curve, the system is allowed to come in and out 

of resonance freely; for the bottom (dashed) curve, the system has only a 20% chance of returning 
to a bound state after leaving. The two dotted curves shows the expected asymptotic behavior of 
the system, i.e., a survival probability of the form Abound oc 
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the orbital elements of extrasolar planets, where e ~ 0.1, the probability of staying in a 3:1 mean 
motion resonance is about 5 times smaller than for a 2:1 resonance. 

2.7. Model Summary 

The basic analytical model developed here contains three physical variables that describe the 
effects of turbulent perturbations on mean motion resonance. The first variable is the natural oscil- 
lation frequency ujq of the resonance (§2.1), modeled here as a pendulum; this quantity represents 
the oscillation frequency of the libration angle in the planetary system. The time scale associated 
with this frequency is typically ~100 times longer than the orbital time scale of the planets, and 
depends on the system properties in a well-known manner (see §2.1 and MD99). The second vari- 
able is the time scale At ~ Stt/O, (or its dimensionless counterpart At w Svrwo/fi) over which the 
turbulent perturbations are re-established to provide an independent realization of the torques. 
The final variable is the amplitude of the perturbations. Since the torques and their corresponding 
changes in angular momentum can be either positive or negative, and hence the mean vanishes, the 
perturbation amplitude can be characterized by the root mean square jrms = ((AJ)^/ J^)^/^; this 
quantity also has a dimensionless counterpart ?7rms given by equations (fT2]) and (fT3l) . Numerical 
simulations of MHD turbulence indicate that the dimensionless amplitude is expected to be of order 
^rms ~ 0.005 (including the reduction factor due to gap clearing). 

For the simplest version of this model, where systems can freely random walk back into a 
resonant state, the three variables (wq, At, r/mis) determine the long term evolution of the ensemble. 
In particular, the expectation value of pendulum energy increases linearly with time according 
to equation ()20p . The probability of the system staying (bound) in resonance decreases as the 
square root of time according to equation ()23|) . Note that both of these results depend only on the 
dimensionless time r = LOot and a dimensionless diffusion constant V = tJj.^J^/At. In other words, 
the action of the turbulent fluctuations can be described by a diffusion constant that determines 
how the variables of the system - considered here as a nonlinear oscillator - random walk through 
phase space. 

In more realistic versions of the model, where systems cannot easily return to a resonant 
state, the long term behavior results in a lower probability of remaining bound. As a result, the 
simple result described above provides an upper limit on the expected number of resonant states 
as a function of time. In the opposite limit in which systems that leave resonance can never 
return, the number of bound states decays exponentially, with varying decay rates, as indicated 
by equations (j26p and (|30p . For the case in which systems can re-enter resonance after leaving, 
but with reduced probability, the number of bound states decreases as t~^/^ as before, but with 
reduced normalization, as shown by equation (j33p . 



3. NUMERICAL SIMULATIONS 



In this section, we consider an ensemble of numerical integrations motivated by an observed 
extrasolar planetary system. For this exploration, we adopt the outer planets c (P ~ 60 d) and b 
(P ~ 30 d) of the Gliese 876 planetary system (Marcy et al. 2001). These planets are observed to lie 
deep in the 2:1 mean motion resonance, with the angles 9i and 02 librating with small amplitudes 
(Laughlin & Chambers 2001, Rivera et al. 2005). Indeed, Gliese 876 b and c represent, by far, the 
most convincing case of an extrasolar planetary system in mean motion resonance. Among the 25 
multiple-planet systems that have been discovered to date, 2:1 mean motion resonances have also 
been tentatively identified for HD 82943 b and c (Gozdziewski Sz Maciejewski 2001, Mayor et al. 
2004, Lee et al. 2006), HD 128311 b and c (Vogt et al. 2005), and HD 73526 b and c (Tinney et 
al. 2005). In each of these three cases, however, the libration widths of the resonant arguments are 
very uncertain, and for HD 128311 and HD 73526, the best dynamical fit to the radial velocity data 
shows only a single argument in libration. Following the discovery of 55 Cancri c and d (Marcy 
et al. 2002), it was suggested by several authors (e.g., Ji et al. 2003) that 55 Cancri b and c are 
possibly participating in a 3:1 mean motion resonance. Self-consistent 6-body fits to the 55 Cancri 
data set published by Fischer et al. (2007), however, show no evidence that 55 Cancri b and c are 
in 3:1 resonance0 

In this set of numerical experiments, we start an ensemble of planetary systems with the obser- 
vationally determined orbital elements of GJ876, so that the system begins in a 2:1 mean motion 
resonance, and then integrate the system forward in time. The integrations include additional 
stochastic fluctuations, which would be present if the circumstellar disk that formed the planets 
were still present. We then monitor how easily the system is knocked out of its resonant configu- 
ration. The ensemble of numerical integrations reported here are thus analogous to those done in 
the previous section with the stochastic pendulum. In this case, however, we integrate the full 18 
phase space variables (6 orbital elements for each planet and 6 more for the star) instead of only 
one variable (93) for the pendulum model. For this system, integrating the motion of the star is 
important because its mass is relatively small, only 0.32 Mq, whereas the masses of the planets 
are 0.79 and 2.53 mj, and hence the system is highly interactive (Laughlin &; Chambers 2001). 
As a result, the GJ876 system is as far from the idealized (one variable) model of the stochastic 
pendulum as any planetary system observed to date; as a result, any agreement found between the 
simple model and the full integrations can be considered robust. 

The starting orbital elements are taken to be those observed (Rivera et al. 2005) at the present 
day; specifically, the epoch for the initial conditions is JD 2449679.6316, the date for which the first 
published Lick Observatory radial velocity was made. The systems are then integrated into the 
future for 2 x 10^ years. Since the planets in this particular system have relatively short periods 
(60.830 days and 30.455 days for the outer and inner planets, respectively), the number of orbits is 



^For an up-to-date x^-ranked list of self-consistent fits to the 55 Cnc data sets see 
links therein. 
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large - about 2.4 x 10^ orbits of the inner planet. The integrations are performed with a Bulirsch- 
Stoer (B-S) integration scheme (e.g., Press et al. 1992) to obtain high accuracy with reasonable 
computational speed. In the absence of turbulent forcing, the system can be integrated forward in 
time for millions of years and is found to stay in resonance, where we use the following definitions 
for the libration angles: 

01 = Al — 2X2 + 2{-cui — vj-i) and ^2 = Ai — 2A2 + (roi — ^2) , (35) 

where the Aj are the mean longitudes and the Wj are the longitudes of periapse of the two planets. 
Note these angles are defined to be a linear combination of the standard ones (e.g., Lee & Peale 
2002), and are chosen because they provide a cleaner separation of the apsidal libration and the 
mean motions (J. Chambers, private communication). We measure the degree of resonance by 
monitoring the maximum amplitude of these libration angles over a given interval of time. For 
this purpose, we use a fixed monitoring time of 100 years, which is much longer than the period of 
the oscillations of the resonance, about 3.5 years, so that many oscillation periods are included in 
the determination of the maximum. This maximum libration angle is then tracked as a function 
of time. For the sake of definiteness, we consider the systems to be bound (in resonance) when 
the maximum angle remains less than Q^^ = 90 degrees; when the maximum angle exceeds 0^, = 
90 degrees, we consider the resonance to be compromised. Note that the exact number of bound 
(resonant) states as a function of time depends on the choice of the angular threshold value 0=,,; 
however, the general trends (shown below) are qualitatively similar for any choice of threshold angle 
in the range 90 < 6^* < 175. 

The effects of turbulence are introduced by adding small velocity perturbations at regular time 
intervals. For the sake of definiteness, we take the forcing time interval to be 1/3 yr, which is about 
four times the period of the inner planet, and is consistent with the time required for the turbulent 
torques to be independent (LSA) and with the model equations of §2. In terms of the simple 
pendulum model of the system, variations of the stochastic forcing times could be scaled out of the 
problem (Appendix A) and would effectively add width to the distribution of forcing strengths. In 
these experiments, we take the size of the velocity perturbations to have the form Aw = va^, where 
the amplitude va = 0.0032 km/s, and where ^ is a uniformly distributed random variable on the 
interval — 1 < C ^ 1- (This amplitude value was chosen to be a round number, namely 2 x 10~^, 
in code units, where G = 1, mass is measured in Mq, and time is measured in years.) Both the x 
and ij components of the velocity are perturbed (independently) in this manner, but the vertical 
velocity component is left unchanged. Including both components of the velocity perturbation, we 
find that the fractional change in speed, and hence angular momentum, has amplitude jrms ~ 
consistent with the benchmark value of equation Q. 

The results from one ensemble of simulations are shown in Figures H] and [5l Figure [4] shows 
the time evolution of the maximum libration angle (as defined above) for five different trials, i.e., 
five different realizations of the turbulent fluctuations. Note that the systems can re-enter bound 
states - defined here to be maximum libration angles less then 6^ = 90 degrees - after leaving. 
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Fig. 4. — Time evolution of the resonance angles for a collection of planetary systems based on the 
observed system GJ876. The systems are subjected to turbulent forcing as described in the text. 
The five curves show the first resonance angle as a function of time (given here in years) for the 
five different realizations of the turbulent fluctuations. 
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Fig. 5. — Time evolution of the fraction of bound states, as measured by the maximum excursions 
of the resonance angles, for a collection of planetary systems based on the observed system GJ876. 
The systems are subjected to turbulent forcing as described in the text. The solid and dot-dashed 
curves show the fraction of systems that remain in resonance as a function of time (in years) for 
the first and second resonance angles, respectively. The dotted curve shows the fraction of bound 
states for the first resonance angle that would result if systems that leave resonance could never 
return to a bound state. For reference, the dashed curve shows the power-law behavior oc t~^/'^ 
expected for a simple stochastic pendulum. Note that the survival probability Pb ^ 0.01 if we 
extrapolate this result out to time scales of ~1 Myr. 
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Nonetheless, the results show a clear trend for systems to preferentially leave resonance, rather 
than return, so the number of bound states decreases rapidly with time. 

We have performed an ensemble of 900 integrations of the system, where each numerical 
experiment uses a different realization of the stochastic perturbations due to turbulence (with this 
sample size, root-A^ fluctuations are ~ 1/30 ~ 0.03). Figure [5] shows the fraction of bound states 
as a function of time. The solid and dot-dashed curves show the fraction of systems that remain in 
resonance as a function of time (in years) for the first and second resonance angles, respectively. The 
dotted curves shows the corresponding fraction of bound states that would remain if the systems 
that leave resonance could never return to a bound state. These results clearly indicate that the 
actual fraction of bound states is systematically larger than the fraction of bound states that would 
result if systems were never allowed to return to a resonant state after leaving (compare the dotted 
and solid curves in Figure [5|). In other words, planetary systems can — sometimes — random 
walk back into a resonant configuration after leaving. This behavior is consistent with that found 
for the evolution of individual systems (shown in Fig. S]). Finally, the dashed curve shows the 
power-law behavior "Pbound for reference. Note that the time evolution of the fraction of 

bound states falls well below the prediction of a pure power-law with Pbound ~ The long 

term behavior of the number of bound states seems to follow a steeper power-law, but does reach a 
full exponential decay during the time window studied (which would be expected if systems never 
return to resonance). 

As an order of magnitude check on the time scales, note that the fractional change in velocity 
per stochastic kick has amplitude Av ~ W~'^Vorb, where Vorb is the orbital speed of the inner planet. 
As outlined in §2, however, the energy associated with the potential well of the resonance is much 
smaller than that of orbit. If we denote the change in velocity required to leave resonance as Vres, 
then Vres ~ O.Olforfe, which in turn implies Av ~ O.Olfres. If Av accumulates as a random walk, 
the system thus requires ~ 10^ stochastic forcing kicks to account for enough energy to compromise 
resonance, and this number translates into a time scale of ~ 3000 yr. Indeed, individual systems 
(see Fig. H} and the ensemble (see Fig. [5]) evolve on a comparable time scale. 

4. DISCUSSION AND CONCLUSION 

The main finding of this paper is that mean motion resonance in planetary systems is readily 
compromised through the action of turbulent fluctuations in circumstellar disks. If MRI operates 
and the accompanying turbulent torques are common during the epoch of planetary migration, 
then planetary systems in mean motion resonances are predicted to be rare. Specifically, even for 
the most favorable case in which the system can freely random walk in and out of resonance, only a 
small percentage of the solar systems that produce pairs of planets in a resonant configuration will 
remain in resonance over typical disk lifetimes of ~ 1 Myr if turbulence is present. On the other 
hand, if planetary systems in mean motion resonance are found to be common, then this analysis 
implies that turbulence, and hence MRI, would have a severely limited duty cycle. 
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We have addressed this problem using both direct numerical integrations of sample solar sys- 
tems and through the construction of model equations that represent an ensemble of stochastically 
driven pendulums. This latter approach allows for the derivation of a number of analytic results 
that elucidate the basic physics of the problem — solar systems being driven from their resonant 
states by turbulent forcing. For this class of systems, the energy of the ensemble of stochastic 
oscillators increases linearly with time (see Fig. [H eq. [20], and Appendix B), the distribution of 
energy for the ensemble spreads according to the solution of equation (fT9]) , and the probability of 
a planetary system remaining in resonance decreases as (see Fig. [2]). This law holds both for 

the case in which systems can freely re-enter resonance, and for the case in which they re-enter res- 
onance with a given (fixed) probability, although the normalization is lower for the latter case (eq. 
[33]). If systems that leave resonance can never return to a bound state, then the time evolution 
of the fraction of bound states is a decaying exponential (§2.5). 

The results from our numerical integrations indicate that the fully interactive system (with 
18 phase space variables) behaves in a qualitatively similar manner to the stochastic pendulum 
(with one variable). In particular, the number of bound states (systems in resonance) decreases 
with time, where the evolution of the ensemble of systems is diffusive, and where the diffusion 
constant depends on the square of the fluctuation amplitude and inversely on the time scale of the 
turbulent forcing. The resulting function that describes the decrease in the number of resonant 
states is intermediate between the power-law behavior of the simplest pendulum system and 
the decaying exponential form that applies for systems that can never re-enter resonance after 
becoming unbound. The numerical simulations show that systems can indeed random walk back 
into a resonant state after leaving, consistent with the finding that the number of bound states 
does not decay exponentially. On the other hand, the return to resonance is relatively rare, in 
particular more unlikely than for the (one variable) stochastic pendulum, a result that is consistent 
with the fraction of bound states falling below the law. This latter power-law behavior — 

the solution to the the long term behavior of a stochastic pendulum - thus provides an upper limit 
on the probability of resonant states as a function of time. 

The most important astronomical implication of this work is a quantitative estimate of the 
rarity of solar systems surviving in a resonant configuration. Our results can be summarized by 
the following expression, which represents our working estimate for the fraction of surviving bound 
(resonant) states: 



where pret is the probability of returning to resonance, r/rms is set by the amplitude of the fluctua- 
tions, and N^rh is the number of orbits over which the turbulence is active. In the second equality we 
have deflned the dimensionless quantity C = (32/7r)^/^/9i.et/??rms- For typical turbulent amplitudes, 
r/rms ~ 0.005, and hence C 640prct- Given that returning to resonance requires somewhat special 
conditions, we expect pret to be small, i.e., p^et *S 1- This expectation is borne out in our numerical 
simulations (§3), which show that "Pbound decreases rapidly with time. As a result, we expect that 




(36) 
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the constant C will be of order ~ 10 (perhaps as large as ~ 50) so that T^bound ~ lO/A'orb ■ 
If turbulence is present for a substantial fraction of the disk lifetime, then A'^oj-b = 10^ — 10*^ and 
T'bound = 0.001 - 0.01, i.e., mean motion resonances will be rare. 

In addition to the prediction that resonant states are rare, this work has other astronomical 
implications: For the particular case of the Gliese 876 system, if no eccentricity damping is assumed 
for the epoch of planet formation and migration, then the current eccentricities and libration widths 
of the system suggest that the system migrated only a small distance (8% of the initial semi-major 
axes) while in resonance (Lee & Peale 2002). This result is consistent with any turbulence in the 
disk of the GJ876 system having had little time or ability to provide substantial perturbations on 
the resonant arguments. On the other hand, the observed HD 128311, HD 82943 and HD73526 
planetary systems are all consistent with having received significant exposure to turbulent per- 
turbations, but not enough to have destroyed their librations completely. Once a large sample of 
multiple-planet systems has been assembled, we should gain further insight from the distribution of 
observed resonant widths, in addition to the observed fraction of planets in mean motion resonance. 
In particular, if planets can random walk back into the resonance after leaving, as indicated by both 
our numerical integrations and model equations, the observational sample should show a definite 
preference for systems with large libration widths. This trend is exactly what we observe for the 
tenuous resonances in the three systems other than Gliese 876. 

This analysis also outlines the basic physical mechanism by which turbulence affects mean 
motion resonance. We have modeled the process through a pendulum equation (which represents 
the resonance) with the addition of stochastic forcing (which represents the turbulence). This 
formulation applies specifically to the class of systems for which the turbulent fluctuations provide 
forcing kicks that are fully independent, i.e., where both the amplitude of the fluctuations and 
the time intervals are not correlated. In this case, the results are largely independent of the 
distributions, and depend primarily on the expectation values of the amplitudes and time intervals 
of the turbulent forcing. These latter quantities set the value of the effective diffusion constant 
(for an analogous problem, see Fatuzzo & Adams 2002, which considers ambipolar diffusion with 
turbulent fluctuations in molecular cloud cores). 

As an intuitive check on the results presented herein, consider the following heuristic argument. 
In a stochastic pendulum, the system accumulates angular momentum as a random walk, so the 
change in angular momentum after A'step steps is given by A^step(AJ)^, where (AJ)^ is the typical 
amplitude of the angular momentum fluctuations. In order of magnitude, the amount of angular 
momentum in the resonant system (when bound) is given by ( Jorb^^o/^)^- The combination of these 
two expressions indicates that the criterion for compromising resonance can be written in the form 
-^step > [(AJ)fc/Jorb]~^('^o/^)^- In Order of magnitude, (AJ)fc/Jorb ~ 10^^ and ujq/VL ~ 10^^, and 
hence the required number of steps is of order A'^step ~ 10^. Indeed, when the number of stochastic 
steps (measured here by the number of orbits A'orb) exceeds this limit, the probability of remaining 
in a bound state decreases (see eqs. [23] and [36]). 
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Notice also that the effective diffusion constant D oc (AJ)^/(At)fe. In order to set the proper 
value of the constant D, one must choose the time interval to be long enough so that the fluctu- 
ations are independent and short enough that the changes in angular momentum (AJ) increase 
linearly with time (where we now suppress the index k). Suppose, for example, one were to 
choose a longer time interval to sample the fluctuations, say (At) ^ (At)o, where (Ai)o is the 
time interval for optimal independent sampling. The change in angular momentum over one 
sampling interval would then increase as a random walk (even within the time interval), i.e., 
(A J) = (A J)o[(At)/(At)o]^^^. The effective value of the diffusion constant then scales accord- 
ing to D oc (AJ)2/(At) = (AJ)g[(At)/(At)o]/(At) = (AJ)g/(At)o. We thus recover the correct 
value, provided that the changes in angular momentum are computed properly. In this paper, we 
have taken the time interval At to be four times the orbital period of the inner planet (twice the 
period of the outer planet), consistent with the calculations of LSA (see also Nelson 2005). The 
above considerations imply that the diffusion constant D is inversely proportional to At, and the 
fraction of bound states scales as Z)^^/^, so that our quoted results are only weakly dependent on 
any uncertainty in the specification of At. 

Although our general conclusion is robust - namely, that turbulence easily compromises mean 
motion resonances - a good deal of additional work should be done. The parameter space available 
for these planetary systems is enormous and further exploration should be carried out for different 
solar system architectures. In particular, the degree to which the full numerical integration agree 
with the simple model equations should be considered, including a determination of the solar system 
properties rcqiiircd for consistency. This treatment assumes that the planets have already formed, 
and have already migrated to their (nearly) final radial locations, and then considers the effects 
of turbulence on mean motion resonances. In actuality, planets can be subjected to turbulent 
fiuctuations as they form, and as they migrate inward and become locked into resonance. These 
processes of planet formation, migration, resonance locking, and turbulent forcing should thus 
be studied concurrently. In addition, the differences between the results of our full numerical 
integrations (which include all 18 phase space variables of the system) and those of the idealized 
stochastic pendulum (which has only one degree of freedom) poses a mathematically interesting 
issue for further study. 
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A. RESCALING FOR VARIABLE TIME INTERVALS 

In this Appendix, we show that if the stochastic forcing impulses occur at irregular time 
intervals, then the time variable can be rescaled to make the stochastic forcing kicks occur at 
regular (periodic) intervals. As a compensating result, however, the distribution of the forcing 
strengths becomes wider and the natural oscillation frequency of the system takes on a distribution 
of values. 

We begin by writing the stochastic pendulum equation in the form 

^ + [cjg + qkQiXkt)] sin = , (Al) 

where Q is a periodic function and qk determines the forcing strength for the kth time interval; the 
variable X/. has unit mean, but allows for different time intervals for the stochastic forcing terms. 
Notice that we generally use periodic Dirac-delta functions to specify the forcing time, but this 
argument can be generalized to include any periodic function Q. If we re-scale the time variable 
for each cycle according to 

t^Xkt, (A2) 
then the stochastic pendulum equation takes the form 

+ H + QkQit)] sin(^ = , (A3) 

where the time variable appearing in this equation is the rescaled one. If we then re-scale the 
remaining parameters (a;o,9fc) according to the relations 

too ^ uJo/Xk = ujk and Qk ^ Qk/ = Qk , (A4) 

the stochastic pendulum equation has the same form as that for the case of perfectly periodic forcing 
terms. This result is thus the analog of Theorem 1 of Adams & Bloch (2008) for the stochastic 
Hill's equation. Notice that in this case the forcing strength acquires a new (generally wider) 
distribution to accommodate the variation in forcing periods, and, the natural oscillation frequency 
lvq becomes a random variable. 



B. ALTERNATE DERIVATION OF THE LINEAR TIME DEPENDENCE FOR 

THE MEAN ENERGY 

In this Appendix, we present an alternate derivation of the result that the energy of a stochas- 
tically driven pendulum tends to increase linearly with time (in the limit of long times). This result 
is thus consistent with equation ()20p in the text. 

We begin by considering the pendulum equation with a periodic delta-function forcing term, 
as described in the text (§2). The angle ip itself must be continuous at the moment of forcing, but 
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the derivative of i-p obeys a jump condition that can be written in the form 



dip 



dip 



(Bl) 



where we have denoted the value of the angle at the kth cycle boundary as (pk- Since the angle 
ip must be a continuous function across the boundary, the potential energy is also a continuous 
function, and hence the total energy experiences discrete jumps of the form 



AE 



-qk sm ipk 



dip 



(B2) 



Over long times, the first term in the above equation averages to zero, whereas the second term 
averages to (g|)/4. The mean energy thus grows as 



4At' 



(B3) 



where At is the time interval of the stochastic forcing. Not only does this argument reproduce the 
linear growth of the mean energy, as found using the Fokker-Planck analysis in the text, but it also 
results in the same coefficient if we identify the effective diffusion constant D = (g^)/(Ar). 
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